Chaotic Systems

Teng-Jui Lin

Content from UW AMATH 301, Beginning Scientific Computing. Although not covered in Spring 2020, the topic is presented in previous years: (1) Application of Runge-Kutta to Lorenz Equation, (2) Vectorized Time-step Integrators.

• Chaotic system
  • Lorenz system
    • One initial condition
    • Two close but distinct initial conditions
    • A cube of initial conditions
  • Visialization techniques
    • Interactive plotting
    • Animations

Lorenz system

The Lorenz system is a system of nonlinear ODEs developed by Edward Lorenz in 1963 to model atmospheric convection. The system

$$ \begin{aligned} \dot{x} &= \sigma (y-x) \\ \dot{y} &= x(\rho-z) - y \\ \dot{z} &= xy - \beta z \end{aligned} $$

with parameters $\sigma = 10, \rho = 28, \beta = 8/3$ can lead to chaotic behavior, where the behavior is unpredictable and sensitive to initial condition and error accumulation.

One initial condition

Problem Statement. Solve the Lorenz system with the initial condition of $x(0) = -8, y(0) = 8, z(0) = 27$ for $t \in [0, 100]$ with $\Delta t = 0.01$ using the above parameters.

• Generate a static plot of all the trajectory in the time interval.
• Generate a rotating animation of all the trajectory in the time interval.
• Generate an animation of the trajectory over time.

Note on animation: for reproducible results, download ffmpeg and add to path variable, then change animation.ffmpeg_path in the plot setting rcParams as your path.

Solving the system and static plot

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy
from scipy import integrate

# Lorenz's chaotic parameters
sigma = 10
rho = 28
beta = 8/3

# initial conditions
initial_val = np.array([-8, 8, 27])

# time array
t_initial = 0
t_final = 100
dt = 0.01
t = np.arange(t_initial, t_final+dt/2, dt)
t_len = len(t)

# define ode system
dxdt = lambda x, y, z : sigma*(y - x)
dydt = lambda x, y, z : x*(rho - z) - y
dzdt = lambda x, y, z : x*y - beta*z
ode_syst = lambda t, v : np.array([dxdt(*v), dydt(*v), dzdt(*v)])

# solve ode system
ode_soln = scipy.integrate.solve_ivp(ode_syst, [t_initial, t_final], initial_val, t_eval=t).y

# plot settings
%config InlineBackend.figure_format = 'retina'
%matplotlib inline

plt.rcParams.update({
    'font.family': 'Arial',  # Times New Roman, Calibri
    'font.weight': 'normal',
    'mathtext.fontset': 'cm',
    'font.size': 18,
    
    'lines.linewidth': 2,
    
    'axes.linewidth': 2,
    'axes.spines.top': False,
    'axes.spines.right': False,
    'axes.titleweight': 'bold',
    'axes.titlesize': 18,
    'axes.labelweight': 'bold',
    
    'xtick.major.size': 8,
    'xtick.major.width': 2,
    'ytick.major.size': 8,
    'ytick.major.width': 2,
    
    'figure.dpi': 80,
    'savefig.dpi': 300,
    
    'legend.framealpha': 1, 
    'legend.edgecolor': 'black',
    'legend.fancybox': False,
    'legend.fontsize': 14,
    
    'animation.html': 'html5',
    'animation.ffmpeg_path': r'C:\\Softwares\\ffmpeg\\bin\\ffmpeg.exe',
})

# plot 3d visualization
plt.rcParams.update({'xtick.labelsize': 18,})
fig, ax = plt.subplots(figsize=(5, 5), subplot_kw={"projection": "3d"})
ax.plot(*ode_soln, '.', color='black', alpha=0.5)

# plot settings
ax.set_xlabel('$x$', fontsize=18, labelpad=20)
ax.set_ylabel('$y$', fontsize=18, labelpad=20)
ax.set_zlabel('$z$', fontsize=18, labelpad=20)
ax.set_title(f'Lorenz\'s Equation \n $\sigma={sigma}, \\beta={beta :.1f}, \\rho={rho}$', fontsize=18, y=1)
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.view_init(30, -45)  # viewing angle
ax.dist = 13  # zoom out to fit all labels in

Rotating animation of all trajectory

# plot static portion
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"projection": "3d"})
ax.plot(*ode_soln, '.', color='black', alpha=0.5)
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.set_axis_off()

def draw_frame(n):
    '''
    Commands to update parameters.
    Here, the viewing angle changes each frame.
    '''
    ax.view_init(elev=0, azim=360*n/(frame_num))
    return fig,

# create animation of given time length
# note here we fit all the data points into the given animation time
from matplotlib import animation
anim_time = 4  # s
fps = 60
frame_num = int(fps * anim_time)
anim = animation.FuncAnimation(fig, draw_frame, frames=frame_num, interval=1000/fps, blit=True)

plt.close()  # disable showing initial frame

# convert animation to video (time-limiting step)
from IPython.display import HTML
HTML(anim.to_html5_video() + '<style>video{width: 400px !important; height: auto;}</style>')

Animation of trajectory over time

# interactive plot
%matplotlib qt

# plot static portion
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"projection": "3d"})
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.set_axis_off()

# plot empty framework
points, = ax.plot([], [], [], '.', color='black', alpha=0.5)
title = ax.set_title('')

# update changes each frame
for n in range(frame_num):
    time_points = round(t_len/frame_num)
    points.set_data(ode_soln[:2, :time_points*n+time_points])
    points.set_3d_properties(ode_soln[2, :time_points*n+time_points])
    title.set_text(f't = {t[n] :.3f}')
    plt.pause(0.0001)

# back to static plot and animations
%matplotlib inline

# plot static portion
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"projection": "3d"})
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.set_axis_off()

# plot empty framework
points, = ax.plot([], [], [], '.', color='black', alpha=0.5)
title = ax.set_title('')

def draw_frame(n):
    '''
    Commands to update parameters.
    Here, the 3d data points changes each frame.
    '''
    num_points = round(t_len/frame_num)
    points.set_data(ode_soln[:2, :num_points*n+num_points])
    points.set_3d_properties(ode_soln[2, :num_points*n+num_points])
    title.set_text(f't = {t[n] :.3f}')
    return fig,

# create animation of given time length
# note here we fit all the data points into the given animation time
from matplotlib import animation
anim_time = 4  # s
fps = 60
frame_num = int(fps * anim_time)
anim = animation.FuncAnimation(fig, draw_frame, frames=frame_num, interval=1000/fps, blit=True)

plt.close()  # disable showing initial frame

# convert animation to video (time-limiting step)
from IPython.display import HTML
HTML(anim.to_html5_video() + '<style>video{width: 400px !important; height: auto;}</style>')

Two initial conditions

Because Lorenz system with the above parameter is a chaotic system, it is sensitive to changes in initial conditions and error accumulation. We show this using the trajectory of two close but distinct initial conditions.

Problem Statement. Solve the Lorenz system with the two initial conditions of

• $x(0) = -8, y(0) = 8, z(0) = 27$
• $x(0) = -8.01, y(0) = 8, z(0) = 27$

for $t \in [0, 100]$ with $\Delta t = 0.01$ using the above parameters.

• Generate a static plot of all the trajectory in the time interval.
• Generate an rotating animation of all the trajectory in the time interval.
• Generate an interactive plot of the trajectory over time. Save the the interactive plot as an animation.

Solving the systems and static plot

# Lorenz's chaotic parameters
sigma = 10
rho = 28
beta = 8/3

# initial conditions
initial_val1 = np.array([-8, 8, 27])
initial_val2 = np.array([-8.01, 8, 27])

# time array
t_initial = 0
t_final = 100
dt = 0.01
t = np.arange(t_initial, t_final+dt/2, dt)
t_len = len(t)

# define ode system
dxdt = lambda x, y, z : sigma*(y - x)
dydt = lambda x, y, z : x*(rho - z) - y
dzdt = lambda x, y, z : x*y - beta*z
ode_syst = lambda t, v : np.array([dxdt(*v), dydt(*v), dzdt(*v)])

# solve ode system
ode_soln1 = scipy.integrate.solve_ivp(ode_syst, [t_initial, t_final], initial_val1, t_eval=t).y
ode_soln2 = scipy.integrate.solve_ivp(ode_syst, [t_initial, t_final], initial_val2, t_eval=t).y

# plot 3d visualization
plt.rcParams.update({'xtick.labelsize': 18,})
fig, ax = plt.subplots(figsize=(5, 5), subplot_kw={"projection": "3d"})
ax.plot(*ode_soln1, '.', alpha=0.7)
ax.plot(*ode_soln2, '.', alpha=0.1)

# plot settings
ax.set_xlabel('$x$', fontsize=18, labelpad=20)
ax.set_ylabel('$y$', fontsize=18, labelpad=20)
ax.set_zlabel('$z$', fontsize=18, labelpad=20)
ax.set_title(f'Lorenz\'s Equation \n $\sigma={sigma}, \\beta={beta :.1f}, \\rho={rho}$', fontsize=18, y=1)
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.view_init(30, -45)  # viewing angle
ax.dist = 13  # zoom out to fit all labels in

Rotating animation of all trajectory

# plot static portion
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"projection": "3d"})
ax.plot(*ode_soln1, '.', alpha=0.7)
ax.plot(*ode_soln2, '.', alpha=0.1)
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.set_axis_off()

def draw_frame(n):
    '''
    Commands to update parameters.
    Here, the viewing angle changes each frame.
    '''
    ax.view_init(elev=0, azim=360*n/(frame_num))
    return fig,

# create animation of given time length
# note here we fit all the data points into the given animation time
from matplotlib import animation
anim_time = 4  # s
fps = 60
frame_num = int(fps * anim_time)
anim = animation.FuncAnimation(fig, draw_frame, frames=frame_num, interval=1000/fps, blit=True)

plt.close()  # disable showing initial frame

# convert animation to video (time-limiting step)
from IPython.display import HTML
HTML(anim.to_html5_video() + '<style>video{width: 400px !important; height: auto;}</style>')

Interactive plot and animation of trajectory over time

# interactive plot
%matplotlib qt

# plot static portion
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"projection": "3d"})
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.set_axis_off()

# plot empty framework
points1, = ax.plot([], [], [], '.', alpha=0.7)
points2, = ax.plot([], [], [], '.', alpha=0.1)
title = ax.set_title('')

# update changes each frame
for n in range(frame_num):
    time_points = round(t_len/frame_num)
    points1.set_data(ode_soln1[:2, :time_points*n+time_points])
    points1.set_3d_properties(ode_soln1[2, :time_points*n+time_points])
    points2.set_data(ode_soln2[:2, :time_points*n+time_points])
    points2.set_3d_properties(ode_soln2[2, :time_points*n+time_points])
    title.set_text(f't = {t[n] :.3f}')
    plt.pause(0.00001)

# back to static plot and animations
%matplotlib inline

# plot static portion
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"projection": "3d"})
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.set_axis_off()

# plot empty framework
points1, = ax.plot([], [], [], '.', alpha=0.7)
points2, = ax.plot([], [], [], '.', alpha=0.1)
title = ax.set_title('')

def draw_frame(n):
    '''
    Commands to update parameters.
    Here, the 3d data points and title changes each frame.
    '''
    time_points = round(t_len/frame_num)
    points1.set_data(ode_soln1[:2, :time_points*n+time_points])
    points1.set_3d_properties(ode_soln1[2, :time_points*n+time_points])
    points2.set_data(ode_soln2[:2, :time_points*n+time_points])
    points2.set_3d_properties(ode_soln2[2, :time_points*n+time_points])
    title.set_text(f't = {t[n] :.3f}')
    return fig,

# create animation of given time length
# note here we fit all the data points into the given animation time
from matplotlib import animation
anim_time = 4  # s
fps = 60
frame_num = int(fps * anim_time)
anim = animation.FuncAnimation(fig, draw_frame, frames=frame_num, interval=1000/fps, blit=True)

plt.close()  # disable showing initial frame

# convert animation to video (time-limiting step)
from IPython.display import HTML
HTML(anim.to_html5_video() + '<style>video{width: 400px !important; height: auto;}</style>')

The above three plots/animations demonstrates that even with very small change in initial conditions, the chaotic nature of the Lorenz system gives vastly different trajectories.

A cube of initial conditions

Problem Statement. Solve the Lorenz system with cubes of initial conditions of

• a cube of points with dimension of 40 with intervals $\Delta x = \Delta y = \Delta z = 2$ centered around $C = (-8, 8, 27)$
  • for $t \in [0, 2]$ with $\Delta t = 1/300$
• a cube of points with dimension of 20 with intervals $\Delta x = \Delta y = \Delta z = 1$ centered around $C = (0, 0, 0)$
  • for $t \in [0, 5]$ with $\Delta t = 1/200$
• a cube of points with dimension of 2 with intervals $\Delta x = \Delta y = \Delta z = 0.1$ centered around $C = (0, 0, 0)$
  • for $t \in [0, 25]$ with $\Delta t = 1/50$
• a cube of points with dimension of 2 with intervals $\Delta x = \Delta y = \Delta z = 0.1$ centered around $C = (1, -1, 5)$
  • for $t \in [0, 15]$ with $\Delta t = 1/100$
• a cube of points with dimension of 2 with intervals $\Delta x = \Delta y = \Delta z = 0.1$ centered around $C = (3, 3, 3)$
  • for $t \in [0, 15]$ with $\Delta t = 1/70$

using the above parameters. Generate an interactive plot of the trajectory over time. Save the the interactive plot as an animation.

Solution. Here, we demonstrate the last set of initial conditions.

Solving the systems and static plot

# Lorenz's chaotic parameters
sigma = 10
beta = 8/3
rho = 28

# initial conditions defining params
dimension = 2
interval = 0.1
center = np.array([3, 3, 3])

# calculate cube of initial conditions
half_dim = dimension / 2
xvec = np.arange(-half_dim, half_dim+interval/2, interval)
yvec = np.arange(-half_dim, half_dim+interval/2, interval)
zvec = np.arange(-half_dim, half_dim+interval/2, interval)
x_len = len(xvec)
y_len = len(yvec)
z_len = len(zvec)
initial_val_vec = np.array(np.meshgrid(xvec+center[0], yvec+center[1], zvec+center[2]))

# time array
t_initial = 0
t_final = 15
dt = 1/70
t = np.arange(t_initial, t_final+dt/2, dt)
t_len = len(t)
t_len

1051

# define ode system
dxdt = lambda x, y, z : sigma*(y - x)
dydt = lambda x, y, z : x*(rho - z) - y
dzdt = lambda x, y, z : x*y - beta*z

def ode_syst(t, v):
    '''System of ODEs for Lorenz system.'''
    # convert 1d input from solve_ivp() back to proper dimension
    v = v.reshape(3, x_len, y_len, z_len)
    return np.array([dxdt(*v), dydt(*v), dzdt(*v)])

# solve ode system
# note: solve_ivp() only accepts 1d array as initial condition, 
# so ode_syst need to covert the initial condition vector back
ode_soln_raw = scipy.integrate.solve_ivp(ode_syst, 
                                     [t_initial, t_final], 
                                     initial_val_vec.reshape(-1), 
                                     t_eval=t, vectorized=True).y
# reshape the soln dimension to proper form
ode_soln = ode_soln_raw.reshape(3, x_len*y_len*z_len, len(t))

# plot 3d visualization
plt.rcParams.update({'xtick.labelsize': 18,})
fig, ax = plt.subplots(figsize=(5, 5), subplot_kw={"projection": "3d"})
ax.plot(*ode_soln[:, :, 0], '.', color='black', alpha=0.1)

# plot setting
ax.set_xlabel('$x$', fontsize=18, labelpad=20)
ax.set_ylabel('$y$', fontsize=18, labelpad=20)
ax.set_zlabel('$z$', fontsize=18, labelpad=20)
ax.set_title(f'Initial conditions', fontsize=18, y=1)
ax.set_xlim(-30, 30)
ax.set_ylim(-30, 30)
ax.set_zlim(0, 50)
ax.view_init(45, -45)  # viewing angle
ax.dist = 13  # zoom out to fit all labels in

Interactive plot and animation of trajectory over time

# animation parameters
# note that here we assume one frame per data dt
fps = 60
anim_time = t_len / fps  # s
frame_num = t_len

# interactive plot
%matplotlib qt

# plot static portion
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"projection": "3d"})
ax.set_xlim(-30, 30)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.set_axis_off()

# plot empty framework
points, = ax.plot([], [], [], '.', color='black', alpha=0.1)
title = ax.set_title('')

# update changes each frame
for n in range(frame_num):
    time_points = round(t_len/frame_num)
    points.set_data(ode_soln[:2, :, time_points*n:time_points*n+time_points].reshape(2, -1))
    points.set_3d_properties(ode_soln[2, :, time_points*n:time_points*n+time_points].reshape(1, -1))
    title.set_text(f't = {t[n] :.3f}')
    plt.pause(0.00001)

C:\Softwares\Anaconda\Anaconda\lib\site-packages\mpl_toolkits\mplot3d\proj3d.py:189: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  return np.array([xs, ys, zs, np.ones_like(xs)])

# back to static plot and animations
%matplotlib inline

# plot static portion
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"projection": "3d"})
ax.set_xlim(-30, 30)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 50)
ax.set_axis_off()

# plot empty framework
points, = ax.plot([], [], [], '.', color='black', alpha=0.1)
title = ax.set_title('')


def draw_frame(n):
    '''
    Commands to update parameters.
    Here, the 3d data points and title changes each frame.
    '''
    time_points = round(t_len/frame_num)
    points.set_data(ode_soln[:2, :, time_points*n:time_points*n+time_points].reshape(2, -1))
    points.set_3d_properties(ode_soln[2, :, time_points*n:time_points*n+time_points].reshape(1, -1))
    title.set_text(f't = {t[n] :.3f}')
    return fig,

# create animation of the data
# note that here we assume one frame per data dt
# the expected animation time is t_len/60
from matplotlib import animation
anim = animation.FuncAnimation(fig, draw_frame, frames=frame_num, interval=1000/fps, blit=True)

plt.close()

# convert animation to video (time-limiting step)
from IPython.display import HTML
HTML(anim.to_html5_video() + '<style>video{width: 400px !important; height: auto;}</style>')

C:\Softwares\Anaconda\Anaconda\lib\site-packages\mpl_toolkits\mplot3d\proj3d.py:189: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  return np.array([xs, ys, zs, np.ones_like(xs)])